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Abstract 

This paper deals with the numerical modeling of wave propagation in porous 
media described by Biot's theory. The viscous efforts between the fluid and 
the elastic skeleton are assumed to be a linear function of the relative velocity, 
which is valid in the low-frequency range. The coexistence of propagating 
fast compressional wave and shear wave, and of a diffusive slow compres- 
sional wave, makes numerical modeling tricky. To avoid restrictions on the 
time step, the Biot's system is splitted into two parts: the propagative part 
is discretized by a fourth-order ADER scheme, while the diffusive part is 
solved analytically. Near the material interfaces, a space-time mesh refine- 
ment is implemented to capture the small spatial scales related to the slow 
compressional wave. The jump conditions along the interfaces are discretized 
by an immersed interface method. Numerical experiments and comparisons 
with exact solutions confirm the accuracy of the numerical modeling. The 
efficiency of the approach is illustrated by simulations of multiple scattering. 
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1. Introduction 



The propagation of waves in porous media has crucial imphcations in 
many areas, such as the characterization of industrial foams, spongious bones 
and petroleum rocks. The most widely used model describing the propaga- 
tion of mechanical waves in a saturated porous medium was proposed by 
Biot in 1956. A major achievement in Biot's theory was the prediction of a 
second (slow) compressional wave, besides the (fast) compressional wave and 
the shear wave classically propagated in elastic media. 

Two regimes are distinguished, depending on the frequency of the waves. 
At frequencies smaller than a critical frequency fc, the fluid flow inside the 
pores is of Poiseuille type, and the viscous efforts between the fluid and 
the solid depend linearly on the relative velocity. In this case, the slow 
compressional wave is almost static and highly attenuated Q. An adequate 
modeling of this diffusive mode remains a major challenge in real applications. 
At frequencies greater than f^, inertial effects begin to dominate the shear 
forces, resulting in an ideal fiowprofile except in the viscous boundary layer, 
and the slow wave propagates j5i_32|. Experimental observations of the slow 
wave in the low-frequency range 36| and in the high-frequency range 10| have 
confirmed Biot's theory. In the current study, we focus on the low-frequency 
range. 

Until the 1990 's, Biot's equations were mainly studied in the harmonic 
regime. Various time-domain methods have been proposed since, based on 
finite differences 
methods 



14, 46, 45 , finite elements 47 , discontinuous Galerkin 



38|, boundary elements [411, pseudospectral methods [8| and spec 



tral element methods 
can be found in j^. 



34j |. A recent review of computational poroelasticity 
Nevertheless, none of the methods proposed in the 



low-frequency range give a complete answer to the following difficulties: 

• the viscous effects greatly influence numerical stability, imposing a re- 
strictive time step. In some physically relevant cases, computations 
cannot be carried out in a reasonable time; 



the wavelength of the slow compressional wave is much smaller than 
that of the other waves. Consequently, one faces the following alter- 
native: either a coarse grid well-suited to the fast wave is chosen, and 
the slow wave is badly discretized; either a fine mesh is used, and the 
computational cost increases dramatically; 
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• maximum computational efficiency is obtained on a Cartesian grid; in 
counterpart, the interfaces are coarsely discretized, which yields spuri- 
ous solutions. Alternatively, unstructured meshes adapted to the inter- 
faces provide accurate description of geometries and jump conditions; 
however, the computational effort greatly increases, due to the cost of 
the mesh generation and to the CFL condition of stability. 

The aim of the present study is to develop an efficient numerical strategy 
to remove these drawbacks. A time-splitting is used along with a fourth-order 



ADER scheme 42| to integrate Biot's equations. A flux- conserving space- 



time mesh refinement |3| is implemented around the interfaces to capture 



the slow compressional wave. Lastly, an immersed interface method |26l . |27 
is developed to provide a subcell resolution of the interfaces and to accu- 
rately enforce the jump conditions between the different porous media. As 
illustrated by the simulations, the combination of these numerical methods 
highlights the importance of an accurate modeling of the slow wave. 



is 



This article, which generalizes a previous one- dimensional work [12 
organized as follows. Biot's model is briefly recalled in section 2. The nu- 
merical methods are described in section 3. Section 4 presents numerical 
experiments and comparisons with exact solutions. In section 5, conclusions 
are drawn and future perspectives are suggested. 



2. Physical modeling 

2.1. Biot's model 

Biot's model describes the propagation of mechanical waves in a porous 
medium consisting of a solid matrix saturated with fluid circulating freely 
through the pores [3, 0, 0, [oj. It is assumed that 

• the wavelengths are large compared with the diameter of the pores; 

• the amplitudes of perturbations are small; 

• the elastic and isotropic matrix is fully saturated by a single fluid phase. 

This model relies on 10 physical parameters: the density pf and the dynamic 
viscosity rj of the fluid; the density ps and the shear modulus p of the elastic 
skeleton; the porosity < < 1, the tortuosity a > 1, the absolute perme- 
ability K, the Lame coefficient Xf and the two Biot's coefficients /3 and m of 
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the saturated matrix. The unknowns are the elastic and acoustic displace- 
ments Us and Uf, the elastic stress tensor a, and the acoustic pressure p. In 
one hand, the constitutive laws are: 

a = (\f tie — l3m^)I + 2ne, 

(1) 

p = m {—(3 tr e + C,) , 
where I is the identity, ^ is the rate of fluid change, and e is the strain tensor 

e = -V. (0(u^-u,)), e = l (Vu, + Vuf). (2) 

The symmetry of a in ([T]) implies compatibility conditions between spatial 



derivatives of e, leading to the Beltrami-Michell equation 39|, Il3 



oxoy OX OX OX oy oy Oy 

a Ao Ao + 2/i 

t'o = ~TT1 : 7? t'l = — ; -, tl2 



(3) 



(4) 



4(Ao + /i)' 4(Ao + ^)' 2(Ao + /i)' 

where Aq = A/ — m is the Lame coefficient of the dry matrix. 
On the other hand, the conservation of momentum yields 

d\s dw 1] 

where Vg = ^ = (wsi, f ^2)^ is the elastic velocity, and w = ^ (uj — u^) = 
(wi, W2Y is the filtration velocity. To be valid, the second equation of (jl]) 
requires that the spectrum of the waves lies mainly in the low-frequency 
range, involving frequencies lower than 

Zti an pf 

If / > fc, niore sophisticated models are required jsl, |29[. In practice, the 
viscosity of the fluid is always non-zero; nevertheless, considering 77 = can 
be relevant for two reasons: 



if / ^ /c, the viscous forces are smaller than the inertial forces |l7l . |34 
and can be neglected to a flrst approximation; 
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• the exact solutions of poro-elastodynamic equations are computed more 
accurately if the saturating fluid is inviscid, which is attractive to val- 
idate the numerical methods. 

2.2. Evolution equations 

A velocity-stress formulation is followed: from ([1]) and (jl]), we obtain the 
system of PDEs 

dan dai2\ pj dp pf r] 



' dvsi 


Pw 


dt 


X 


dVs2 


Pw 


dt 


X 


dwi 
dt 


+^ 

X 


dW2 

dt 
< „ 


+^ 

X 



dx dy J X dx x ^ 
d(7i2 da22\ Pf dp Pf T] 

H T: TT- — UJ2, 



dx dy J X dy x 



dx dy J X dx x ^ 
^cri2 da22\ , p dp p 7] 

\- I H — = W2, 



dx dy J X dy x « 

dan , „ ,dvsi „ dwi dvs2 n dw2 , 

— — - (A/ + 2/i) — pm— Xf — Pm—— = U^^, 

dt dx dx dy dy 

doi2 f dVs2 dvsi\ , 



(6) 



dt \ dx dy 

d(J22 . dvsi n dwi .dVs2 o dw2 , 

-ITT- - (^/ + 2^) -t: P^^r~ = J' 

dt dx dx dy dy 

dp [ ndvsi , dwi ndvs2 , dw2\ „ 
— + m /3 — — + — — + {3 — — + — — = /p, 
dt \ dx dx dy dy J 

where /^,,, f^,^, f^^^ and fp are force densities, p^ = ^pf^ p = 0p^+(i-0) p^, 
and X = P Pw — p} > ^- Setting 



0"22 ' 



U = {vsi, Vs2, wi,W2, au, ai2, (J22, pY , 

F = (0, 0, 0, 0, /ctu, /cri2) /ct22 5 fp) 1 

equations (E]) are written as a first-order non-homogeneous linear system 



(7) 



^U-haAu + bAu=-SU + F, (8) 



where A, B and S are 8x8 real matrices detailled in [Appendix A 
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The energy of poroelastic waves can be deduced from (jS]). Without any 
source terms (F = 0) and setting C e = cr + it is proven in |16i] that 

E{t) = l[ {p^rs'+Ce■.e)dS+l- [ {p^w^ + -p^) dS+ [ pf^r,.w dS (9) 

is an energy that satisfies 

I' "^-^-dS. (10) 



d t J ^2 K 

Consequently, E is conserved when the viscous effects are neglected [r] = 0) 
and is a decreasing function otherwise. 

2.3. Heterogeneous media 




n 

Figure 1: Interface F between two poroelastic media JIq and fii 

The physical parameters defined in section [2Jl are piecewise constant and 
can be discontinuous across interfaces. In what follows, we will focus on two 
domains Qq and Qi, which are separated by a stationary interface F described 
by a parametric equation (a;(r), ?/(r)) (figured]). At any point P on F, the 
unit tangential vector t and the unit normal vector n are 

x\y] , n = ^=^= fy', -x') . (11) 



The derivatives x' = and y' = ^ are assumed to be continuous everywhere 
along F, and to be differentiable as many times as required. 

The evolution equations ([6]) must be completed by a set of jump con- 
ditions. The simple case of perfect bonding and perfect hydraulic contact 



along F is considered here, modeled by the jump conditions [18 



[v,] = 0, [w.n] = 0, [a.n]=0, [p] = 0. (12) 
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Enforcing these conditions is one of the main objective of the immersed 
interface method presented in section 13.31 




Figure 2: Phase velocities (a) and attenuations (b) of the solutions to Biot's model cor- 
responding to the porous medium Oq in table [21 pf: fast compressional wave; ps: slow 
compressional wave; s: shear wave. In (a), the horizontal dotted lines refer to the eigen- 
values Cpf, Cps and Cg of A and B. 

The eigenvalues of A and B in ([H]) are reah ±Cp/, ±Cps, ±Cs, and 
(multiphcity 2), where Cpj > max(cs, Cps) > 0. If ?7 7^ 0, the spectral radius 
R{S) = ^ ^ can be very large and then the system ([8]) is stiff. 

A plane wave d e*^"^ is injected in ([6]), where k = /ce and d are the 
wavevector and the polarization, respectively; r is the position, a; = 2 tt / is 
the angular frequency and / is the frequency. If d is coUinear with k, the 
dispersion relation of compressional waves is obtained: 

Ak^ + B{u) k^ + C{u) =0, 

A = Km [Xf + 2 fi — (3'^ m) , C{co) = x f^^^ — iv P^^^ (13) 
B{u) = -K{{Xf + 2 fi) p^ + m{p-2l3 pf))u'^ + iri{Xf + 2^)u, 

where the roots ±/Cp/ and ±kps satisfy < ^e{kpf} < ^e{kps}. If d is 
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orthogonal with k, the dispersion relation of the shear wave is obtained: 



k = 



AC-B^^ 1/2 



A = u'^ {p + (p Pf{a - 2)) - iu <f'^ 



B = —uj^(ppf {a — 1) + iuj (p'^ —, 



(14) 



C = CO (j) Of a — iuj (j) 

where the roots are ±ks, ^e{ks} > 0. Based on (fT3|) and (IT^ . the phase 
velocities c = u/^e{k} and the attenuations a = 53m {/c} of each wave are 
defined. In the remainder of this article, the subscripts pf, ps and s denote 
fast compressional, slow compressional and shear waves, respectively. 

The phase velocities Cpf{f), Cps{f) and Cs(/) are monotonically increasing 
functions, tending asymptotically towards the eigenvalues Cpf, Cpg and Cg- If 
rj = 0, the three waves are non-dispersive and non-dissipative, and the energy 
of poroelastic waves iQ is conserved. If r/ 7^ 0, the fast compressional wave 
and the shear wave are weakly dispersive and dissipative. The slow compres- 
sional wave, however, is highly modified by the viscosity of the saturating 
fluid. If / ^ fc, then Cps{f) <^^, and the slow compressional wave tends 
towards a static diffusive mode [lO|. At greater frequencies, Cps is larger but 



the attenuation increases. These properties are summarized in figure [2l 

In the low-frequency range, the direct contribution of the slow compres- 
sional wave to the overall wave propagation processes is therefore negligible 
when considering an homogeneous medium. However, the influence of the 
slow wave becomes crucial in heterogeneous media Q. The slow compres- 
sional wave, generated during the interaction between the propagative waves 
and the scatterers, remains localized around the interfaces. Consequently, 
this slow wave has a major influence on the balance equations at the in- 
terfaces, modifying crucially the behavior of fast compressional and shear 
diffracted waves. An accurate computation of the slow wave is therefore 
necessary, as shown in the numerical tests of section |H 



3. Numerical modeling 

3.1. Numerical scheme 

To integrate the system (|8]), a uniform grid is introduced; the spatial mesh 
sizes are Ax, Ay, and the time step is A t. A straightforward discretization 
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of ([H]) by an explicit time scheme leads to the following stability condition: 



eAx 2 \ 

Cpf 



^'^™"l^'i?(S)j' '''' 



where G is obtained by a Von-Neumann analysis of stability when S = 0. In 
(fTSjl . the bound induced by the spectral radius of S can be very restrictive. In 
sandstone saturated with bitumen, for example, the maximal CFL number 
is roughly CpfAt/Ax ^ 10^^^ ^ G, which is intractable for computations. 
We follow here a more efficient strategy based on second-order Strang's 



splitting 23|, solving alternatively the hyperbolic system 



^U + aAu + bAu = 0. (16) 

and then the diffusive system with a source term 

|-U = -SU + F. (17) 

ot 

The linear system ( fT6l) is solved by applying any scheme for hyperbolic sys- 
tems, giving U"j^^^^. In the numerical experiments performed in section HJ a 



fourth-order ADER scheme is used 42 , which involves a centered stencil of 



25 nodes. On Cartesian grids, this scheme amounts to a fourth-order Lax- 



Wendroff scheme [28| . It is dispersive of order 4 and dissipative of order 6 



and its stability limit is G = 1 |43l . |25| . Other single-grid schemes can be 
used without any restrictions. 

Since the physical parameters do not vary with time, the diffusive system 
( |T7l) is solved exactly. For simplicity, null force density is taken: F = 0. In 



this case, p and a are unchanged, whereas the velocities become (fc = 1, 2) 

n+1 n+1/2 , Pf -RPt\ n+1/2 

p \ ) ^ (18) 

„+l -HPT n+1/2 

where T depends on the time step (see section [331). The splitting (IT6l)- (|T711 
along with exact integration (ITSl) recovers the optimal condition of stability: 

CpfAt/Ax < G. 

Since the matrices A and B do not commute with S, the theoretical 
order of convergence falls from 4 to 2 when the viscosity is non-negligible. 
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Using a fourth-order accurate scheme such as ADER 4 is nevertheless advan- 
tageous, compared with a second-order scheme such as Lax-Wendroff: the 
stabihty hmit is improved, and numerical artifacts (dispersion, attenuation, 
anisotropy) are greatly reduced. 



In [38[ , the authors notice that the first-order splitting does not lead to a 
correct representation of the slow mode at low fre quen cies. Nevertheless, the 
numerous one-dimensional examples provided in [l2| demonstrate that the 
second-order splitting accurately represents the static mode when a sufficient 
number of discretization points per wavelength is used. This can be obtained 
by using a local space-time refinement presented in the following section. 

3.2. Mesh refinement 

The slow wave has much smaller spatial scales of evolution than the wave- 
length of the other waves. A very fine grid is therefore required to account for 
its evolution. Since the use of a fine uniform grid on the whole computational 
domain is out of reach, grid refinement provides a good alternative. In ad- 
dition, the slow wave remains localized near the interfaces (section 12.41) . and 
hence grid refinement is necessary only around these places. Lastly, even if 
the slow wave propagates {rj = 0) the property Cpf 3> Cps is usually satisfied: 
consequently, a fine mesh near the interface is still useful to perform accurate 
extrapolations, as required by the immersed interface method (section l3.3p . 

We adopt here a space-time mesh refinement approach based on fiux 
conservation [sf, which is more naturally coupled to the fiux-conserving 
scheme developed to solve f lT6|) . The refined zones are rectangular Cartesian 
patches with mesh sizes Ax / q, Ay / q, where the integer q is the refinement 
factor. To reduce the cpu time and to limit the numerical dispersion on 



the coarse grid, a local time step At / q is used [22|, |40| . When one time 
step is done on the coarse grid, q time substeps are done on the refined 
zone. The extrapolated values required to couple coarse and fine grids are 
obtained by linear interpolation in space and time on the numerical values at 
the surrounding nodes j25|. In the case of the Lax-Wendroff scheme applied 
to the scalar advection equation, the stability of the coupling is proven in ^ 
whatever q. 

The additional cost induced by mesh refinement can become prohibitive, 
both concerning the memory requirements and the computational time, be- 
cause of the q substeps inside one time step. The value of q must be estimated 
carefully in terms of the physical parameters. For this purpose, the wave- 
lengths Xpfifo) = Cp/(/o)//o and Aps(/o) = Cps(/o)//o are deduced from the 
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dispersion analysis, where /o is the central frequency of the source. The num- 
ber of fine grid nodes per wavelength of the slow compressional wave and the 
number of coarse grid nodes per wavelength of the fast compressional wave 
must then be roughly equal: 

V(/o) ^ Kfifo) ^ ^ ^ Cpfifo) 
Ax / q Ax Cpsifo)' 

3.3. Immersed interface method 

The discretization of the interfaces requires special care. A straightfor- 
ward stair-step discretization of the interfaces introduces a first-order geo- 
metrical error and yields spurious numerical diffractions. In addition, the 
jump conditions ( JT2l) are not enforced numerically if no special treatment is 
applied. Lastly, the smoothness requirements to solve f|T6l) are not satisfied, 
decreasing the convergence rate of the ADER scheme. 

To remove these drawbacks while maintaining the efficiency of Cartesian 
grid methods, we adapt an immersed interface method previously developed 



in acoustics and elastodynamics [35|, |26|, |27|, |25[ . At the irregular points where 



the ADER's stencil crosses the interface F, the scheme will use modified values 
of the solution, instead of the usual numerical values. The modified values 
are extrapolations, based on the local geometry of F and on r successive 
derivatives of the jump conditions ( !T2ll . The parameter r is discussed at the 
end of this section. 

Let us consider a point M{xi, yj) G Qi and its orthogonal projection P 
onto F (figure [3]). The algorithm to build the modified value at M is divided 
into four steps. 



Step 1: high-order interface conditions. 

On the side [k = 0, 1), the boundary values of the spatial derivatives of 
U up to the r-th order are put in a vector U)^ with = 4 (r + 1) (r + 2) 
components: 

UI,= lim fu^..., ^ , ^\ U^,..., #-U^V, (20) 

where / = 0, r and m = 0, /. Following this formalism, the zero-th 
order jump conditions (IT^ are written 

C?U? = C°U°, (21) 
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where the matrices of the jump conditions C° depend on the local geometry 
of r: 



1 1 




















0\ 
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1 / 



(22) 



The jump condition (pT]) is differentiated with respect to time t, and then the 
time derivatives are replaced by spatial derivatives thanks to the conservation 
law (ITB]) . For example, we obtain 



■cSA„AuS 



d 



(23) 



where Aq and Bq are the matrices in VLq. The jump condition (]2Ti) is also 
differentiated in terms of r. Taking advantage of the chain-rule, we obtain 
e.g. 



d ^0 



UO + Co(.'^UO + ,^Uo). (24) 



From ([21]), ([23D and ([24]), we build matrices such that C^U} 
which provides first-order jump conditions. By iterating this process r times, 
r-th order interface conditions are obtained 



C\\J\ = Cl\]l (25) 

where CJJ are Uc x n„ matrices {k = 0, 1), and = 3 (r + 1) (r + 2). The 
computation of matrices is a tedious task when r > 2, that can be greatly 
simplified using computer algebra tools. 



Step 2: high-order Beltrami-Michell equations. 

The equation ([3]) is satisfied anywhere in a poroelastic medium. Under suffi- 
cient smoothness requirements, it can be differentiated with respect to x and 
as many times as required: 



d'j (Ti2 



d x^ ' ^ dy 



i+l 



00 



d^a22 , . d^p 



dx^ * 9 



dx^ * 9 



dx^ * (9 

ail . n 9^ '^22 , ^ P 



(26) 
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where j > 2 and i = 0, ■ ■ ■ , j — 2. The equations are also satisfied 
along r. They can be used therefore to reduce the number of independent 
components in U^. For this purpose, we define the vectors such that 



(27) 



where G)^ are Uy x [riy — 71^) matrices, and = r (r — l)/2 if r > 2, rif, = 
otherwise. Based on fl2U]) and fl2^ . an algorithm to compute the non-zero 



components of G^. is proposed in [Appendix B 



Step 3: high-order boundary values. 

Based on ( |25|) and (!27|) . the vectors of independent boundary values satisfy 



(28) 



where 



are Uc x {riy 



rib) matrices. Since the system ( 125]) is 
underdetermined, the solution is not unique, and hence it can be written 



v[ = ((SD-' SSlKs.) 



A" 



(29) 



where (S^)~^ is the least-squares pseudo-inverse of S^, Ksj is the matrix 
filled with the kernel of S^, and A^ is a set of — Uc — Ub Lagrange multi- 
pliers that represent the coordinates of V[ onto the kernel. A singular value 
decomposition of is used to build (S^^)"^ and the kernel Kgr [37 . 



Step 4: construction of modified values. 

Let be the matrix of r-th order 2D Taylor expansions 



n 



1, 



1 



m] 



iVj - VpY 



T ! 



I. 



(30) 
The 



where Is is the 8x8 identity matrix, I = 0,..., r and m = 0,..., /. 
modified value at (x/, yj) is a smooth extension of the solution on the other 
side of r (figure [3]), and it writes: 



U 



i,j 



-rrr j jr 



(31) 



The vector Vq in fl3T]) remains to be estimated in terms of the boundary 
conditions and of the numerical values at surrounding grid points. For this 
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I 



Figure 3: Irregular point M{xi, yj) e fli and its orthogonal projection P onto F. The 
grid nodes used to compute j are inside the circle with radius d and centered on P; 
they are denoted by +. 



purpose, we consider the disc V centered at P with a radius d, that contains 
Nd grid points. At the grid points of P fl f2o, '^-th order Taylor expansion of 
the solution at P gives 



r+l\ 



V5\ ^ (32) 

m,Gai|o) . +0(Ax^+i). 



A"" 



At the grid points of V (iQi, r-th order Taylor expansion of the solution at 
P and the boundary conditions ( l29l) give 

/ vj; \ (33) 
= ni^Gi {{Sir' S5|Ks.) ( \+0{Ax^+'). 



14 



Equations f l32p and fl33p are written in the matrix form 



(U(., Q) 



V 



M 



+ 



(34) 



where M is a convenient 8 N^x [2 — 2 rih — matrix. To ensure that the 
system fl5^ is overdetermined, the radius d of the disc is chosen to satisfy 



e{d, r) 



2 72^ — 2 nf, 



> 1. 



(35) 



Exact values in are replaced by numerical ones, and the Taylor rests are 
removed. The least-squares inverse of M is denoted by M^^ The La granere 
multipliers A^ are accounted in the construction of M, but are not involved 
in the definition of the modified value fl3Tl) . As a consequence, they can be 
removed and the [n^ — rih) x 8 restriction of is defined by 

V 



Lastly, the modified value follows from ( 13T|) and ( 136|) : 

U 



n^jG^M-i (U" 



(36) 



(37) 



Comments and practical details. 



1. A similar algorithm is applied at each irregular point along F. The 
sizes of the matrices involved are summarized in table [TJ Since the 
jump conditions do not vary with time, the evaluation of the matrices 
in f l37|) is done during a preprocessing step. Only small matrix-vector 
products are therefore required at each time step. After optimization 
of the computer codes, this additional cost is made negligible, lower 
than 1% of the time-marching. 

2. The matrix M in depends on the sub cell position of P inside the 
mesh and on the jump conditions at P, involving the local geometry 
and the curvature of F at P. Consequently, all these insights are incor- 
porated in the modified value ( 137|) . and hence in the scheme. 
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u-OjIJ-ui u y 


Size 


riy 


4(r + 1) (r + 2) 




ric 


3 (r + 1) (r + 2) 




rib 


r(r- 1)/ 2 if r > 2, 


else 


CI 






CI 


riv X (n„ - Hb) 






ric X (n^ - rib) 






8 X riy 




1 


{riy - rib) X {n^ - rib] 







{riy - rib) X (n^ - Uc 


- nb) 


M 


8NdX (2 - 2 rib - 


ric) 




(n^ - rib) x8 Nd 





Table 1: Quantities involved in the computation of the modified values (section [3.3p . 



3. The simulations indicate that overestimation of e in (135|) has a crucial 
influence on the stability of the immersed interface method. Various 
strategies can be used to ensure (|35ll . for instance an adaptive choice 
of d depending on the local geometry of F at P. We adopt here a 
simpler strategy, based on a constant radius d. Taking r = 2, numerical 
experiments have shown that d = 3.2 Ax is a good candidate, while 
d = 4.5 A a; is used when r = 3. In this case, we obtain typically 
Nd ^ 20 and e ^ 4. 

4. The order r plays an important role on the accuracy of the coupling 
between the immersed interface method and a s-th order scheme. If r > 
s, then a s-th order local truncation error is obtained at the irregular 
points. This condition can be slightly relaxed: r = s—1 still ensures a s- 
th order overall accuracy [19|]. As a consequence, a fourth-order ADER 
scheme (s = 4) requires a third-order immersed interface method (r = 
3) to maintain fourth-order convergence. 

5. A GKS analysis of stability has been performed in ID in the case of 



an inviscid saturating fluid [25[. Extending this approach to 2D prob- 
lems with viscous saturating fluids is out of reach. Various numerical 
experiments, however, indicate the stability of the method under the 
usual CFL condition (section [3. ip . if two requirements are satisfied: (i) 
the number of grid nodes used for extrapolations is sufficiently large, 
as stated in point [3l (ii) the Beltrami-Michell equations (1271) are used. 
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3.4- Summary of the algorithm 

The numerical strategy presented in this section couples three numerical 
methods: a finite difference numerical scheme with splitting (section l3.ip . a 
space-time mesh refinement (section [3.2p . and an immersed interface method 
(section l3.3p . To clarify the interactions between these methods, the global 
algorithm is summarized as follows: 

> Preprocessing 

- Detection of irregular grid points 

- Computation of extrapolation matrices in (I37D 

- Initialization of the solution at t = 

- Diffusive step G18I) where T = At / 2 on the coarse grid 

- Diffusive step G18D where T = At/(2q) on the refined grids 

> Time iterations 

• Coarse grid: 

- Computation of modified values (\37\) if present 

- Solving the propagative step G16D 

- Diffusive step G18D where T = At 

• On each refined grid, q subtime iterations: 

- Space-time interpolations at the grid boundaries 

- Computation of modified values G37D 

- Solving the propagative step G16D 

- Diffusive step G18D where T = At / q 

> End of time iterations 

- Diffusive step G18I) where T = At / 2 on the coarse grid 

- Diffusive step G18I) where T = At/{2q) on the refined grids 
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4. Numerical experiments 

4.I. Configurations 

Five tests are proposed along this section. In Test 1, the convergence 
order of the ADER scheme coupled with the immersed interface method 
is measured. Test 2 illustrates the different kind of waves in homogeneous 
media, and also the influence of the local space-time refinement. Test 3 
investigates the numerical stability of the global algorithm. Diffraction of a 
plane wave by one (Test 4) and four (Test 5) cylindrical scatterers illustrates 
the accuracy and the physical relevance of the proposed numerical methods. 

The physical parameters given in table H] correspond to Cold Lake sand- 
stone and shale saturated with water [l^, respectively. In some experiments, 
an inviscid saturating fluid is artificially considered: t] = Pa.s, the other 
parameters being unchanged. As recalled in section 12. ![ this limit-case has 
physical significance only in the high-frequency range. It is mainly addressed 
here for a numerical purpose. 



Parameters 


^^0 


9.1 


Ps (kg/m^) 


2650 


2211 


jj (Pa) 


2.926 10^ 


3.539 10^ 


Pf (kg/m^) 


1040 


1040 


r] (Pa.s) 


1.510-3 


10-3 





0.335 


0.05 


a 


2 


2 


K (m^) 


10-" 


5. 10-12 


Xf (Pa) 


6.1425 10^ 


4.689 10^ 


(3 


0.9558 


0.0527 


m (Pa) 


6.491 10^ 


9.852 10^ 


Cpf (m/s) 


2384.1 


2350.4 


Cps (m/s) 


758.9 


486.4 


Cs (m/s) 


1229.0 


1290.0 


fc (Hz) 


3844.9 


765.1 



Table 2: Physical parameters of the matrix (ilo) and of the scatterer (fli), corresponding 
to sandstone and shale saturated with water, respectively. 

Once the spatial mesh sizes A x and A y are chosen on the coarse grid, 
the time step follows from the CFL number in Cp/o At/ max (A x, Ay) = 
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0.95 < 1. If ?7 7^ 0, the maximum CFL number induced by fllSp is equal to 
0.5: consequently, the simulations done here with the splitting (fTB]) - (fT7|) are 
twice faster than with unsplitted methods. 

The grids are excited by two means: either a plane fast compressional 
wave, either a point source that generates cylindrical waves. Details of the 
excitation method are given in Appendix C In the case of an incident plane 
wave, the exact expression given in Appendix C is also enforced numerically 
on the edges of the computational domain. No special attention is paid to 



simulate outgoing waves, for instance with Perfectly-Matched Layers [3l|, |46 



In all the presented tests, the size of the domain and the number of iterations 
are chosen to avoid the spurious reflections of diffracted waves with the outer 
frontiers. 



4-2. Test 1: convergence measurements 

In Test 1, we focus on the coupling between the ADER scheme (section 
13.11) and the immersed interface method (section 13.31) . For this purpose, we 
consider a domain [0, 400] m^ cut by a plane interface with slope 60 degrees. 
The saturating fluids are inviscid: exact solutions can be computed very 
accurately without Fourier synthesis, and splitting errors of the scheme are 
avoided. The source is plane (]C.2I) . with parameters: 7 = 10~^, /o = 40 Hz, 
t = 3.3 10~^ s, and 9 = —30 degrees (figure IH-a) . Consequently, the incident 
wave propagates normally to the interface, leading to a 1-D configuration; 
from a numerical point of view, however, the problem is fully bidimensional. 

The computations are done on a uniform grid of N x N points, during 
A^/4 time steps. Comparisons with the exact values of the pressure p are 
done on the subdomain [50, 350] m x [150, 250] m, in order to avoid spurious 
effects induced by the edges of the computational domain (figure HJ-b). The 
measures involve reflected and transmitted fast and slow compressional waves 
generated by the interface (no shear wave is generated in ID); these waves 
are highly sensitive to the discretization of the jump conditions. 

Errors in I2 norm and convergence rates are reported in tableland drawn 
on figure O Various values of r are investigated; r = means that no 
immersed interface method is applied; in this case, first-order accuracy is 
obtained. As stated in section 13.31 fourth-order accuracy is maintained if 
r = 4 — 1 = 3, i.e if third-order extrapolations are used in the immersed 
interface method. In the present test case, r = 2 is sufficient to obtain the 
same level of accuracy on a large range of grid size. Nevertheless, this could 



19 




Figure 4: Test 1. Snapshots of p at the initial instant (a) and at the instant of measure 
(b). The white rectangle denotes the zone where convergence errors are measured. 



be untrue in other contexts, and hence we will always use r = 3 in the 
following simulations. 

4-3. Test 2: mesh refinement 

In the second test, we focus on the coupling between the ADER scheme 
(section 13. ip and the mesh refinement (section 13. 2p . For this purpose, a 
homogeneous medium Qq on a domain [—250, 250] m^ is excited by the force 
density ( ]C.4p . The parameters of the source are: Xg = Us = 0, /o = 40 Hz, 



C = ll^, ro = 2(^, and 7 = 10^. The computational domain is discretized 
on a coarse mesh of 500^ points. Two locally refined areas are added: one 
around the source point = [-25, 25]^, and one at ^2 = [80, 130] x [80, 120]. 
Both grids are refined by a factor g = 5, leading to 255^ points in Qi and 
205 X 255 points in Q2. 

Figure [6] shows snapshots of the fields after 280 iterations. At this time 
t ~ 4.4 / /o, the fast compressional wave has also crossed the grid ^2- In the 
left column, one takes 1] = 0, and consequently all the waves generated by 
the source propagate radially. In the right column, the viscosity is non zero, 
and the slow compressional wave remains localized around the source point. 
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N 


r = 




order 


r = 1 




order 


r = 2 




order 


r = 3 




order 


400 


4.894 lO'J 




- 


6.52710" 






6.10710° 




- 


5.06710° 






800 


1.24710° 




1.973 


1.66710° 




1.961 


1.06510° 




2.520 


8.642 10" 


1 


2.552 


1200 


6.520 10" 


1 


1.599 


6.75810" 


1 


2.242 


2.27310" 


1 


3.809 


1.93610" 


1 


3.690 


1600 


4.770 10- 


1 


1.086 


3.61710" 


1 


2.173 


7.995 10" 


2 


3.632 


6.835 10" 


2 


3.619 


2000 


3.818 10" 


1 


0.998 


2.25410" 


1 


2.119 


3.34410" 


2 


3.906 


2.88810" 


2 


3.861 


2400 


3.128 10- 


1 


1.093 


1.50910" 


1 


2.201 


1.56410" 


2 


4.168 


1.41410" 


2 


3.917 


2800 


2.69610- 


1 


0.964 


1.09710" 


1 


2.069 


9.15710" 


3 


3.473 


7.69710" 


3 


3.945 


3200 


2.390 10" 


1 


0.902 


8.431 10" 


2 


1.971 


5.59810" 


3 


3.685 


4.50410" 


3 


4.013 


3600 


2.11510" 


1 


1.038 


6.57310" 


2 


2.114 


3.61310" 


3 


3.718 


2.84010" 


3 


3.915 


4000 


1.92810" 


1 


0.879 


5.36210" 


2 


1.933 


2.48510" 


3 


3.552 


1.86610" 


3 


3.986 



Table 3: Test 1. Convergence rate in I2 norm. No immersed interface method (r — 0), 
linear (r — 1), quadratic (r — 2) or cubic (r ~ 3) immersed interface method. 

This wave is clearly observed on the pressure field, but is also present in the 
stress fields although not visible on the corresponding panels of figure El 

In order to evaluate the infiuence of the refined reference solution 

is computed using a fine mesh Ax = Ay = 1/5 m on the whole computational 
domain. As observed in figure [7]- (a), (b), the local space-time refinement ap- 
plied in Qi and Q2 does not lead to significant spurious refiections during the 
propagation of the waves. On the contrary and as expected, the refinement 
around the source point leads to a much better resolution of the slow wave 
in the viscous case: see figure [7]-(c). 

In conclusion, mesh refinement coupled with the ADER scheme (with or 
without splitting) accurately represents the behavior of the different poroe- 
lastic waves in an homogeneous medium, even though they present a spatially 
complex structure, as observed in figure [7]-(d). 

4-4- Test 3: stability of the complete algorithm 

Some theoretical results based on GKS theory exist concerning the sta- 
bility of immersed interface methods or local mesh refinement [3, 0, 
0, y, [s^- However, these analyses have been done mainly in the case of 
one-dimensional model problems and basic numerical schemes. The present 
algorithm combines more sophisticated numerical methods in 2D, and hence 
GKS analysis is out of reach. The only reasonable way to confirm the stability 
of the full method is to perform numerical perturbation tests. 
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Figure 5; Test 1. Error measured in I2 norm versus the number of points iV^;, with various 
order r of the immersed interface method. Dotted hne corresponds to 4-th order slope. 

For this purpose, we consider a heterogeneous porous media made up by 
a matrix Qq with a cyhndrical scatterer Qi of radius r = 10 m, centered in 
a domain [— 50, 50]m^. In both Qo and fli, the viscosity of the saturating 
fluids are taken into account (77 7^ 0). All the components of U (jTl) are 
initialized randomly at each grid points, and the source F is set to zero. 
The computational domain [—50, 50] m^ is discretized by a coarse grid. The 
scatterer is included in a grid Qi of size 15 x 15 m^ with a large refinement 
factor q = 9. 

In figure [H the energy of poroelastic waves ([9]) is displayed during 3 10^ 
iterations. Theoretically, this energy should slowly decrease, depending on w 
in (ITUI) . At the beginning of the simulations, a large decrease of E is observed. 
It is logically induced by the random non-smooth initial field, which generates 
a large numerical dissipation. After roughly 1000 time steps, a smooth field 
is reached and the mechanical energy slowly decreases. This confirms the 
stability of the full numerical method, involving the ADER scheme with 
splitting, mesh refinement, and immersed interface method. 

4-5. Test 4: diffraction of a plane wave by a cylinder 

A cylindrical scatterer Qi of radius 40 m is centered at point (0, 0) in 
the computational domain [—200, 200] m^. The source is a plane wave (lC.2p 
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Figure 7: Test 2 in the inviscid case, corresponding to the left column of figure [6l Stress 
cTii along lines ?/ = 30 m (a) and y = 130 m (b). Comparison between the fine grid 
reference solution and the solution obtained with a coarse grid and two refined areas of 
factor q = 5 (c). 3D representation of an (d). 
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Figure 8: Test 3. Time evolution of the computed total energy (|9]). 



initially in medium f2o, with parameters: 7 = 1, /o = 40 Hz, t = —2.09 10~^ 
s, and ^ = degree. The initial conditions are illustrated in figure [9] (a-b). 

In figure [HI the viscosity has been canceled in both media Qq and Qi. In 
such a configuration, the analytical solution of ([7]) can be computed using 
Fourier and Bessel expansions and is used to validate the simulations. The 
diffracted waves propagate with velocities Cpf,Cps and Cg given in table |2j To 
ensure the same number of points per wavelength for all the diffracted waves 
(see section [321), the computational domain is locally refined by inserting the 
cylinder in the grid Qi = [—45,45]^ with a refinement factor q = 5, deduced 
from relation (fT9|) and table [2l The infiuence of this refinement combined 
with the immersed interface method is clearly visible in figure [H 

Without refinement nor interface method (g = 1, r = 0, left column), 
the waves created during the interaction with the scatterer are polluted by 
spurious numerical artifacts. With g = 5 and r = 3 (right column), these 
non-physical perturbations disappear and the refiected-transmitted waves 
are correctly computed. Comparisons with the exact solution presented on 
figure [To] confirm the accuracy of the simulation. Without mesh refinement 
nor immersed interface method, inaccurate results are obtained, especially for 
the refiected fast compressional wave where a shift of about 5 m is observed 



In figure [To] -(c) and (d), the time evolution of the pressure registered 
at the point (—60, 60) is presented from t = 0.028 s, in order to avoid the 



(figure [lOl-(b)). 
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Figure 9: Test 4, inviscid fluids, (a): p at initial time, (b): p along the line y = 0; 
vertical lines denote the frontiers of the scatterer. (c-d): p after 180 iterations, (c): no 
refinement nor interface method {q = 1, r = 0). (d): refinement factor q = 5 and third 
order immersed interface method {r — 3). Dotted square represents the frontiers of the 
refined domain. The cross denotes the location of the receiver. 
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Figure 10: Test 4, inviscid fluids, (a): p on the line y = after 180 iterations, (b): zoom 
around the reflected fast compressional wave, (c): time evolution of p recorded at receiver 
{x = —60,y = 60) from t = 0.028 s to the final instant, (d) zoom on the time evolution 
of scattered waves. 
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incident initial wave. We observe the fast reflected wave, followed by the slow 
wave and a combination of diffracted fast and slow waves. Once again, the 
full strategy captures accurately all the temporal variations of the solution, 
while the basic algorithm gives very poor results, especially concerning the 
waves diffracted by the scatterer (see figure [TOl-(d)). 

The same configuration is now considered by taking into account the vis- 
cosity of the saturating fluids (see table [2]). Based on the dispersion analysis 
performed in section 12. 4[ the values of the phase velocity at /o = 40 Hz give 
a refinement factor g 22 in the grid Qi to satisfy our refinement criterion 
( IT^ . Snapshots of p and an after 180 iterations are given on figure [HI 
showing different structures than in the inviscid case presented in figure Hi 
The diffracted fast compressional and shear waves propagate, while the slow 
compressional wave remains localized around the scatterer. 

The spatial structure of this diffusive wave is clearly visible in figure UTl 
(c),(d). Figure [TTI-(e) focuses on the pressure component of the slow wave 
around the interface. It is observed that the numerical method applied with- 
out mesh refinement nor interface treatment gives inaccurate results about 
this static wave. This is also the case of the propagative reflected compres- 
sional wave shown in figure [11]- (f). The computed fields converge when the 
refinement factor in Qi increases. With refinement factors q larger than 15, 
the results are mostly indistinguishable, and hence they are not represented 
here. With q = 15, the refined grid Qi contains 1365^ discretization points 
and involves 19184 irregular points in the immersed interface algorithm. 

The time evolution of the mechanical energy (|9]) is illustrated in figure 
[T2l The times ti = 0.008 s and t/ — 0.066 s correspond respectively to the 
instants when the plane wave begins to interact with the scatterer and when 
it has crossed the scatterer entirely. In the case of inviscid saturating fluids 
(blue circles on figure [T2|) . the energy is almost conserved, as expected, and 
only a decrease of 0.15% is observed. In the viscous case (blue diamond 
on figure [T2|) . the energy decreases when the plane wave interacts with the 
scatterer, and remains constant otherwise. This behavior is logical, since the 
rate of decrease of energy in f IlOp is governed by w. In homogeneous medium, 
w is extremely small compared with other fields, and hence ^ ~ 0. On 
the other hand, during the interaction between the wave and the scatterer, 
the filtration velocity of the slow wave generated at the interface has an 
important amplitude. The theoretical decrease of energy, obtained by a 
numerical integration of equation f llOl) . is shown using a dotted line in figure 
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(a) (b) 
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Figure 11: Test 4, viscous fluids. Snapshot oi p (a) and an (b) after 180 iterations, with a 
refinement factor g = 15 and a third order interface method (r = 3). Snapshot of p around 
the scatterer (c), (d). Zoom oi p on the line y = around the reflected fast compressional 
wave (e) and around the diffusive slow wave around the right interface (f). 



rC2\ and is very close to the observed decrease of E. This confirms that our 
numerical strategy accurately models the dissipation of mechanical energy. 
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Figure 12: Test 4. Time variation of the computed mechanical energy E{t)/E{0) defined, 
by relation ([9]) for inviscid and viscous case. Numerical integration of relation (ITO|) is also 
represented by the dashed curve. 

4-6. Test 5: multiple scatterers 

In the last test, the ability of the proposed numerical strategy to handle 
complex geometries is illustrated. Four cylindrical scatterers of medium Qi 
are inserted in the matrix f^o- Each cylinder is surrounded by a refined grid. 
The refinement factor is g = 5 in the inviscid case, and g = 10 in the case of 
viscous saturating fiuids. The pressure field is represented in figure [T3i The 
behavior of the fast compressional waves is qualitatively the same in both 
cases, unlike the slow waves: 

• in the inviscid case (figure [131 l^ft column) , the slow waves propagate 
and interact with the other scatterers, which generates new sets of 
refiected-transmitted waves. Accurate computation of these successive 
interactions is obtained thanks to the combination of the local refine- 
ment procedure and the immersed interface method; 

• in the viscous case (figure [131 right column) , the slow waves remain 
localized around the interfaces, and hence they do not participate di- 
rectly to the scattering process. At time t2, they have mostly disap- 
peared due to their physically diffusive behavior. As indicated by test 
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Figure 13: Test 5. Snapshots of p at two different times ti ant t2, corresponding to 140 
and 280 iterations. Left column: inviscid fluids; right column: viscous saturating fluids. 
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4, a fine modeling of tliese waves is nevertheless necessary in order to 
compute accurately the propagative diffracted waves. 



5. Conclusion 

Numerical modeling of 2-D poroelastic waves was addressed numerically 
in the time-domain. The evolution equations were issued from Biot's theory, 
assuming viscous efforts of Poiseuille type, which is valid essentially in the 
low-frequency range. Three numerical tools were combined: a fourth-order 
scheme with time-splitting, a space-time mesh refinement, and an immersed 
interface method. The resulting method provides highly accurate simulations 
of wave propagation in realistic configurations. Numerical experiments have 
indicated that even if the slow waves remain localized around sources or in- 
terfaces, the manner in which they are computed influences the propagative 
waves, that are measured in practical applications. Academic cylindrical 
geometries have been considered here; however, more complex smooth ge- 
ometries can be handled without any restrictions, for example cubic splines. 
Various extensions of the present work are suggested: 

• The numerical methods presented here make it possible to simulate rel- 
evant physical experiments. We have especially in mind the modeling 
of multiple scattering in random media. Based on simulated data, the 
properties of the effective medium equivalent to the disordered medium 



under study can be deduced ll|. This numerical approach compares 
advantageously with the methods usually followed by physicists: real 
experiments are expensive, and analytical methods are restricted to 



very small concentrations of scatterers [30[. The example given in sec- 
tion 14.61 is obviously preliminary: interaction of a plane wave with 
hundreds of scatterers needs to be addressed, which requires the par- 
allelization of the algorithms. 

Poroelastic media in perfectly bonded contact are considered here. 
More realistic conditions can be studied and properly enforced by the 
immersed interface method, for instance sliding and imperfect bond- 



ing 27|, or imperfect hydraulic contact (6|. Comparison of numerical 
simulation with experimental results could help to validate or improve 
the models of contact, which constitutes a current issue in poroelastic- 
ity. Generalizing our approach to the interface between a poroelastic 
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medium and a fluid is another direction of work witli large applications, 
for instance in biomechanics. 



Incorporation of attenuation in the elastic skeleton is required to model 
the real processes of dissipation j8|. For this purpose, memory variables 
need to be introduced in the evolution equations (E])- 

Lastly, the numerical modeling of the transient Biot equations in the 
full range of validity of poroelasticity constitutes a natural extension 
of the present work. At frequencies greater than /c in (|5]), a correction 
of the viscosity proportional with the square root of frequency needs 



to be introduced, as described e.g. by the JKD model [2l|, |15[. In the 



time domain, the new evolution equations involve fractional derivatives 
of order 1/2, whose efficient numerical evaluation is a major challenge 



29, 20, 33 . 
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Appendix A. Matrices involved in system ([5]) 

Based on and ([7]), the matrices in (|8]) write 

/ 





A 



(A/ + 2/i) 


-/3 m 








-yu 





-A/ 


-/3 m 





(3 m 


m 






-pw/x -pf/x \ 

-p^/x 

pf/x p/x 

Pf/x 



B 



-p 



-A; 

(A/ + 2p) 









/3 m 


( ° 





-pfix 

















p/x 












-/3 m 


-/3 m 
m 






-pw/x 














-pw/x 


-pf/x 





PfIx 














PfIx 


p/x 






p/x 



0/ 



(A.l) 



Appendix B. Algorithm for the Beltrami-Michell conditions 

The following algorithm is proposed to compute the non-zero components 
of matrices [k = 0, 1) involved in the Step 2 of the immersed interface 
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method fsee section 1331) : 



a = —1, P = —1, 
for 7 = 0, r, for 5 = 0, 7 
if 5 = then for e = 1, 8 

a = a + l, /3 = /3 + l, G^[a,/3] 
if 7 7^ and 5 7^ and 7 7^ 5 then 
if 7 = 2 then u = 0, t] = 0, 
else if (5 = 1 then z/ = 0, 77 = 1, 
else if 5 = 7 — 1 then = 1, 77 = 0, 
else u = 1, T] = 1, 
for £ = 1, 5 

a = a + l, /3 = /3 + l, G^[a,/3] 



a = a + l, /3 = /3-8 + z/, 
/3 = /3 + 2 - z/, 
/3 = /3 + l, 
(3 = (3 + 12, 
P = P + 2-r], 
(3 = (3 + 1, 

a = a + 1, (3 = (3 — 9 + r], 

a = a + 1, (3 = (3 + 1, 
if 7 7^ and 7 = 5 then for e = 1, 

a = a + 1, (3 = (3 + 1, 



Gl[a,P] 
Gl[a,(3] 
Gl[a,(3] 
Gl[a,(3] 
Gl[a,P] 
Gl[a,P] 
Gl[a,P] 
Gl[a,(3] 



1 

Oo 
01 

■00 
02 

1 

■ 1 



(B.l) 



Gl[a,(3] = l. 



Appendix C. Implementation of sources 

Two sources are considered. The first one involves a plane right-going 
fast compressional wave, whose wavevector k makes an angle 6 with the 
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horizontal x-axis. Its time evolution is 



h{t) 



4 I 

am sm{(3m uJot) a <t < — 



m=l 



(C.l) 



otherwise, 



where = 2™" , = 27r/o; the coefficients are: ai = 1, 02 = —21/32, 
03 = 63/768, 04 = —1/512, ensuring smoothness. The support of the 
incident plane wave lies initially in fig- If ^7 7^ 0, this wave is slightly dispersive 
and its time-domain expression follows from a Fourier synthesis: 
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J 



(C.2) 



where 7 is an amplitude factor, kpf is the wavenumber f[T^ and 

((l-0)p, + p//3(a-l)) a;2_(A^ + 2p-m/32) A;2^ + ,c^0/3^ 



pf{a(5- (j))u^ 



luj 



(C.3) 

When ?7 = 0, fcp/ depends linearly on uj, and consequently Ypf and the vector 
column in flC.2p no more involve cu: a straightforward time-domain expression 
of the incident plane wave can be obtained. 

As a second source, we also implement force densities acting on ai2 in 
(E]). The only non-null component in ([7]) is 



/<7i2 = g{x,y)h{t) 



(C.4) 
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where is a truncated gaussian centered at point {xg, Vs)' 

f 7e"(7) if r^y/{x- x^f + {y - Vsf < tq, 
9{x,y)^< (C.5) 
[ otherwise, 

and /i is a Ricker signal: 
y otherwise. 

(C.6) 

This source generates cyhndrical waves of all types: fast and slow compres- 
sional waves, shear waves. 
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